In-class Exercise 3

Author

Ho Zi Jun

Published

September 9, 2024

Modified

September 9, 2024

3 Network Constrained Spatial Point Patterns Analysis

3.1 Loading R packages

pacman::p_load(sf, spNetwork, tmap, tidyverse)

3.2 Data Import and Preparation

The code chunk below uses st_read() of sf package to import Punggol_St and Punggol_CC geospatial data sets into RStudio as sf data frames.

#Punggol_St is in ESRI Shapefile format

network <- st_read(dsn = "data/geospatial",
                   layer = "Punggol_St")
Reading layer `Punggol_St' from data source 
  `C:\zjho008\ISSS626-GAA\In-class_Ex\In-class_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM

Data has to be in linestring and not multiple linestring. Use st to convert it from multi line to single line.

childcare <- st_read(dsn = "data/geospatial",
                     layer = "Punggol_CC") %>%
  st_zm(drop = TRUE,
        what = "ZM") # to remove z value
Reading layer `Punggol_CC' from data source 
  `C:\zjho008\ISSS626-GAA\In-class_Ex\In-class_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM

We can examine the structure of the output simple data features data tables in R studio.

childcare File is in kml format hence the dimension shown is XYZ (additional dimension). It is important to note the dimension. Upon further inspection under geometry the childcare data has point Z.

Note: for take home exercise under entire data file, 1 folder as rawdata, with another separate folder as data for analysis.

3.3 Visualising the Geospatial Data

Before we jump into the analysis, it is a good practice to visualise the geospatial data. There are at least two ways to visualise the geospatial data. One way is by using plot() of Base R as shown in the code chunk below.

plot(st_geometry(network)) # plotting the road network first, especially when in sf layer
plot(childcare, add = T, col = 'red', pch = 19) # followed by the childcare # since mapped with colours when plotted multiple colours do not appear

# add = T -> T = TRUE the point is plotted twice.

Code chunk result when removing the st_geometry:

plot(network)
plot(childcare, add = T, col = 'red', pch = 19)

Note

network has 3 columns: Link ID St_name Geometry

Removing st_geometry will result in individual columns which are pulled out and plotted individually.

To visualise the data with high cartographic quality and in an interactive manner. The mapping function of tmap package can be used as shown in the code chunk below.

tmap_mode('plot')
tmap mode set to plotting
tm_shape(childcare) + # specifying the layer that is being used
  tm_dots(col = "red") +
  tm_shape(network) + # to use the extent of the map layer
  tm_lines()

tmap_mode('plot')
tmap mode set to plotting

ways to add markers;

https://r-tmap.github.io/tmap/reference/index.html

Specify the shape object: tm_symbols() tm_squares() tm_bubbles() tm_dots() - to keep the size constant when performing zoom functions tm_markers()

Making the plot an interactive layer

tmap_mode('view') #just by switching to 'view' to achieve the interactivity
tmap mode set to interactive viewing
tm_shape(childcare) + # specifying the layer that is being used
  tm_dots(col = "red") +
  tm_shape(network) + # to use the extent of the map layer
  tm_lines()
tmap_mode('plot') # to ensure after the session is ended it will end in the plot mode to reduce resource consumption
tmap mode set to plotting
Note

Childcare & network can be switched on and off accordingly.

3 different data consumptions: 2 layers of ESRI map data (WorldGray Canvas & OpenStreetMap) Topographic layer

mpabox ~ leaflet

While using tmap methods requires a longer code chunk the benefit it brings are the flexibility and customisation that can be done.

3.3.1 Preparing the lixels objects

Before computing NKDE, the Spatial Lines object needs to be cut into lixels with a specified minimal distance. This task can be performed by using lixelize_lines() of spNetwork as shown in the code chunk below.

lixels <- lixelize_lines(network,
                         700,  
                         mindist = 350)

given that it is a road network and in the context of the childcare - so using the reasonable walking distance based on weather and perceived hindrance is about 700 metres based on a study for perceivable walking distance.

mindist is set as half for the minimum walking distance.

2642 segments in the line network. to split into line segment each should be 700 and in the centre the minimum distance should be 350m.

After running the code chunk the segments, the remaining is slightly greater than 350.

if increase to 500 the segment is 2645

if reduce to 150m the segment is still the same 2645.

For take home ex 3, the BMR has to be plotted - a rough gauge of the general distance so we should not use a distance smaller than the point. Calculating the nearest neighbour to find out the nearest neighbour -

based on distances starting on the lowest 25 percentile of accidents along the road segment. Want to acheive a segment that can pick up some accident occurences.

3.3.2 Generating line centre points

Next, lines_center() of spNetwork will be used to generate a SpatialPointsDataFrame (i.e. samples) with line centre points as shown in the code chunk below.

samples <- lines_center(lixels) # sf format

3.4 Visualising the lixel segment

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels) +
  tm_lines() +
  tm_shape(samples) +
  tm_dots(size = 0.01)
tmap_mode('plot')
tmap mode set to plotting

3.4.1 Performing NKDE

We are ready to compute the NKDE by using the code chunk below:

densities <- nkde(network,
                  events = childcare,
                  w = rep(1, nrow(childcare)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300,
                  div = "bw",
                  method = "simple",
                  digits = 1,
                  tol = 1,
                  grid_shape = c(1,1),
                  max_depth = 8,
                  agg = 5,
                  sparse = TRUE,
                  verbose = FALSE)

# avoid gaussian if intensity changes to negative # 3 methods: simple, continous, discontinous

the computed density values (i.e. densities) into samples and lixels objects as density field.

samples$density <- densities
lixels$density <- densities

To append the intensity values into the simple tibular frame or lixel data frame simialr to a left join.

Avoid sorting to avoid changing the sequence.

values attached to the line and point.

Since the svy21 projection system is in metres, the computed density values are very small i.e. 0.0000005. The code chunk below is used to rescale the density values from number of events per metre to number of events per kilometre.

# rescaling

samples$density <- samples$density*1000
lixels$density <- lixels$density*1000

The code below uses appropriate functions of tmap package to prepare interactive and high cartographic quality map visualisation.

tmap_mode(mode = c("view"))
tmap mode set to interactive viewing
tm_shape(lixels) +
  tm_lines(col = "density") +
  tm_shape(childcare) +
  tm_dots()
tmap_mode("plot")
tmap mode set to plotting
kfun_childcare <- kfunctions(network,
                             childcare,
                             start = 0,
                             end = 1000,
                             step = 50,
                             width = 50,
                             nsim = 99, # simulations are starting from zero
                             resolution = 50,
                             verbose = FALSE,
                             conf_int = 0.05)

The output of kfunctions() is a list with the following values:

  • plotkA, a ggplot2 object representing the values of the k-function
  • plotgA, a ggplot2 object representing the values of the g-function
  • valuesA, a DataFrame with the values used to build the plots

Visualising the ggplot2 object of k-function by using the code chunk below.

::: panel-tabset ## K-Function

kfun_childcare$plotk # whether to plot G / K function

2 possible patterns observed

regular pattern below the envelope - showing signs of regularity - childcare centres near to each other e.g. at 200m apart which is showing the signs of regularity

and complete spatial randomness at the upper portion.

3.5 G-Function

kfun_childcare$plotg # whether to plot G / K function

both functions are returned